Generalized Manna sandpile model with height restrictions 
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Sandpile models with conserved number of particles (also called fixed energy sandpiles) may 
undergo phase transitions between active and absorbing states. We generalize the Manna sandpile 
model with fixed number of particles, introducing a parameter — 1 < A < 1 related to the toppling 
of particles from active sites to its first neighbors. In particular, we discuss a model with height 
restrictions, allowing for at most two particles on a site. Sites with double occupancy are active, and 
their particles may be transfered to first neighbor sites, if the height restriction do allow the change. 
For A = each one of the two particles is independently assigned to one of the two first neighbors and 
the original stochastic sandpile model is recovered. For A = 1 exactly one particle will be placed on 
each first neighbor and thus a deterministic (BTW) sandpile model is obtained. When A = — 1 two 
particles are moved to one of the first neighbors, and this implies that the density of active sites is 
conserved in the evolution of the system, and no phase transition is observed. Through simulations 
of the stationary state, we estimate the critical density of particles and the critical exponents as 
functions of A. 

PACS numbers: 05.70.Ln, 02.50.Ga, 64.60. Cn,64.60.Kw 



I. INTRODUCTION 

Problems related to phase transitions between active 
and absorbing states have attracted much interest in re- 
cent years [j. Although these transitions occur away 
from thermodynamic equilibrium, since the presence of 
absorbing states prevents detailed balance to be satis- 
fied in the dynamical evolution, the theoretical frame- 
work developed for equilibrium phase transition may in 
fact be applied to these systems, and concepts such as 
scaling and universality are relevant for non-equilibrium 
phase transitions as well. It is then of interest to identify 
the universality classes in these systems. Many of the 
stochastic models with absorbing states which show a 
phase transition belong to the directed percolation (DP) 
universality class including the much studied contact 
process (CP) 0|. The so called DP conjecture con- 
firmed so far through many examples, states that models 
with a scalar order parameter exhibiting phase transi- 
tions between an active state to a single absorbing state 
and without additional conservation laws should belong 
to this universality class. In models for sandpiles the 
number of particles is conserved and an infinite number 
of absorbing states is present, and thus they are poten- 
tial candidates to belong to a non-DP universality class. 
Careful simulations of the unrestricted Manna sandpile 
model @ and the model with height restrictions 0| lead 
to exponents at variance with the DP value. 

The Manna sandpile model was originally proposed as 
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a stochastic model for SOC, with slow addition of sand 
and avalanches which lead to abrupt loss of sand 0|. The 
scaling behavior in the SOC regime was later recognized 
to be associated to an absorbing state phase transition in 
the corresponding model without addition or loss of par- 
ticles (grains of sand) , the so called fixed energy sandpiles 
(FES) Q. In this model, a (i-dimensional lattice of size 
L d is occupied by N particles. A configuration of the lat- 
tice is specified fixing the number of particles Zi at each 
site i, where Zi may be any non-negative integer num- 
ber. Sites with Zi > 2 are active, and an active site may 
lose two particles to its first neighbors, with a unitary 
toppling rate. The two particles move to randomly and 
independently chosen sites among the first neighbors of 
site i. Any state in which no site has two or more par- 
ticles is an absorbing state, and in the thermodynamic 
limit the number of absorbing states is infinite, as long 
as the density of particle £ = N/ L d is smaller than unity. 
Numerical simulations @ in fact show that this model in 
one dimension undergoes a continuous phase transition 
at a critical density Cc = 0.9488 with critical exponents 
which are different from the DP values. 

In this paper we study a variation of the original model 
where the occupancy numbers are restricted to Zi < 2. 
Each move of a particle is accepted only if this con- 
straint is not violated at the destination site. This ad- 
ditional restriction in the allowed configurations leads to 
simplifications in mean-field calculations and in simula- 
tions. Numerical simulations for the one-dimensional re- 
stricted model lead to a slightly lower critical density 
Cc = 0.92965, while the exponents are compatible with 
the ones found for the unrestricted model jfj. We gen- 
eralize this restricted model by introducing a parameter 
whose value is related to the choice of the destination 
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sites of the two toppling particles. The restricted Manna 
model corresponds to A = and a model where exactly 
one particle is sent to each first neighbor of the active 
site (in one dimension) is recovered for A = 1. We call 
this latter case the BTW model, since it corresponds to 
the original Bak, Tang and Wiesenfeld model |H with 
conserved number of particles. In particular, we are in- 
terested in investigating the phase transition in the model 
as the parameter A is changed. 

In the section m we define the model in more details. 
The numerical simulations which were performed are de- 
scribed in section lllll as well as the results they furnished. 
Conclusions may be found in section IWl 



II. DEFINITION OF THE MODEL 

We consider N particles located on the L sites of a 
one-dimensional lattice in such a way that each site i is 
occupied by = 0, 1 or 2 particles. Sites with two par- 
ticles are active and the control parameter of the model 
is the density of particles £ = N/L. The update of the 
configuration is started by randomly choosing one of the 
Na active site of the lattice. The following updates of 
the configuration may then occur: 

1. With a probability (1 + A)/2 each first neighbor 
receives exactly one of the particles from the active 
site, as may be seen in figure 

[p=(l+X)/2] 



Figure 1: The two particles of the active sites are moved to 
different first neighbor sites. 



If one of the first neighbors is active, the particle 
which was chosen to be sent to this site remains 
on the original site. If both first neighbor sites are 
active, no particle movement occurs. This assures 
that the restriction on the occupancy numbers is 
always satisfied. 

2. With a probability (1 — A)/2, one of the two first 
neighbors is chosen with probability 1/2 and both 
particles are moved to this site. If the destination 
site is already occupied by one particle, only one 
of the particles originating from the active site is 
moved. If the destination site is active, no move- 
ment is done. An example of such a transition is 
shown in figure [3 

One iteration as described above corresponds to a time 
increment of 1/Na- After each iteration, the list of active 
sites is updated. In the restricted model the density of 
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Figure 2: Both particles from the active site are moved to the 
same first neighbor site. 



particles obeys £ < 2, and when £ = 2 an additional ab- 
sorbing state is reached. In the unrestricted model there 
is no such upper limit for £ and no spurious absorbing 
state with unitary density of active sites exists. However, 
since the phase transition in the models occurs at particle 
densities below unity, this details do not bother us here. 
The fact that the restriction on the occupancy numbers 
implies a reduction of only about 2% in the critical den- 
sity of the model reveals that for the unrestricted model 
in the supercritical region close to the phase transition 
the fraction of sites with occupancy larger than 2 should 
be very small. 

The parameter A is restricted to the interval [—1,1]. 
In the upper limit of this interval, the toppling of par- 
ticles is deterministic and the model corresponds to a 
BTW sandpile with conserved number of particles. In 
the lower limit the number of active sites is a conserved 
quantity in the time evolution of the model, with active 
sites diffusing on the lattice. The conservation of the 
number of active sites in this limit is a consequence of 
the occupancy numbers restriction and does not happen 
in the unrestricted model. 



III. SIMULATIONS AND RESULTS 

We realized simulations to find the critical properties 
of the one-dimensional model with height restrictions. 
The initial condition, for a given value of A and C = 
N/L, is a uniform and uncorrelated distribution of the 
N particles on the L sites of the lattice, respecting the 
restrictions. We studied lattice sizes L between 100 and 
2000, performing N r repetitions with times up to t max 
for a certain range of densities and each value of A. In 
our simulations t max was in the interval [4 x 10 4 , 2 x 10 7 ] 
and N r = 2000. Since in finite systems with absorbing 
states the only stationary states are the absorbing states 
themselves, to study the transition we are interested in 
the quasi-stationary states of the model. Usually, the 
simulational determination of this state is hindered by 
the presence of the absorbing states, since considerable 
fluctuations will be found in the quantities estimated in 
the simulations, particularly for particle densities close to 
the critical value. To avoid these drawbacks and increase 
the precision of the estimates of the stationary state we 
used the prescription for simulation of quasi-statio nary 
states proposed recently by de Oliveira and Dickman [lf| . 

In the simulations, the time evolution of the density 
of active sites p a is obtained, one example being shown 
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in figure A stationary value of this density may then 
be estimated for given values of ( and A calculating the 
mean value of the last 2000 points in such simulations. 
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Figure 3: Example of a simulation of the time evolution of 
the density of active sites. In this case L = 100, A = 0, and 
C = 0.92. 



For each value of the lattice size L, realizing simula- 
tions for a range of particle densities (, we obtain curves 
for the order parameter such as the one depicted in figure 
Q] It should be mentioned that for each lattice size the 
particle density may assume only a discrete set of val- 
ues, but it seems reasonable to interpolate between those 
values, assuming that p a varies continuously with £. 



sidering the following scaling relations: 

Pa (( c ,L) ~ L-V"±, 

Pa(C) ~ (C-CcY 3 . 



(1) 

(2) 



Here p a (C) denotes the thermodynamic limit L — > oo of 
p a (C,,L). These scaling relations provide estimates for 
the critical exponents (3 and v±_ , as well as for the critical 
particle density Cc- To estimate the critical particle den- 
sity and the order parameter exponent (3 we proceed as 
follows: for a given value of L we choose a critical value 
which maximizes the correlation in a linear approxi- 
mation for the function In p a — b + (3^ ln[£ — C= ] , with 

C > Q • Once obtaining estimates for (3^ and £c L ' ) for 
different values of L, they may be extrapolated to the 
limit L — > oo. The ratio f3/v± may be estimated directly 
through the relation hi p a = c + [3 / v± InL, using values 
for p a at the estimated critical particle density. 

Using the procedure described above, we obtained es- 
timates for the case A = 0, as may be seen in figured We 
then realized that the estimates for (3 may be improved 
excluding the smaller sizes in the extrapolation to the 
thermodynamic limit. A better result for the ratio j3/vx 
is found if we apply the scaling correction proposed in 
for the simulation of the quasi-stationary state of the 
contact process. It seems that the algorithm proposed in 
this reference for simulations of the quasi-stationary state 
and used in our work implies a scaling correction which is 
characteristic of the model. With these corrections, our 
estimates are very close to the ones in reference |fg. The 
extrapolations with the scaling correction may be seen in 
figure El 
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Figure 4: Estimated values of the order parameter p a as a 
function of the particle density £, obtained for L = 2000 and 
A = 0. 



The estimated values of the order parameter as a func- 
tion of the particle density ( may then be analyzed con- 
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Figure 5: Extrapolation without the scaling correction for 
A = 0. The estimates are: Cc = 0.92996, @/v x = 0.2832 and 
P = 0.42608. 



Repeating this strategy for different values of A, we 
finally may obtain the phase diagram of the model, as 
well as study the values of the estimates for the critical 
exponents as functions of this parameter. In figure 
the phase diagram is shown. Besides the points for the 
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Figure 6: Extrapolations with the scaling correction for A = 0, 
Corrected values for the estimates: f c = 0.92996, /3/v± = 
0.247 and (5 = 0.41275. 



critical line which emerge from the simulations, the result 
of a 2-site mean field approximation for the model may 
also be seen in the figure. Details of these calculations 
will be given elsewhere. The curve which results from the 
mean-field calculations is always below the results of the 
simulations, confirming that when the correlations are 
ignored the size of the active region in the phase diagram 
is overestimated. As A approaches the limiting value -1, 
the results of the simulations get closer to the ones in 
the mean-field approximation, which leads to £ c = 1/2 in 
this limit, even in a one-site approximation. This seems 
to be reasonable, since in this limit the evolution of the 
model is dominated by the diffusion of active sites, and 
this leads to a mean-field like behavior. This argument 
applies only in the limit of A slightly above -1, since as 
already noticed above when A = —1, p a (t) is constant 
and no phase transition is observed. 



ical exponents show a rather strong dependence of the 
parameter A, which may be due to a crossover between 
universality classes in the limits A — > — 1, where mean- 
field behavior is expected, and A — > 1, the conservative 
BTW model. In this latter limit, our simulations indi- 
cate a first-order transition at C = 1, as may be seen in 
figure 0J 




Figure 8: Estimates for the critical exponents P e P/vx as 
functions of the parameter A. 
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Figure 7: Phase diagram of the Manna model with height 
restriction parametrized by A. 



In figure |H1 we notice that the estimates for the crit- 
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Figure 9: Estimates for the order parameter for L = 500, with 
A = 1, indicating a first order transition. 



IV. CONCLUSION 

We studied a generalization of the one-dimensional 
Manna model with height restrictions and conservation 
of the number of particles, with the inclusion of a pa- 
rameter A which is related to the two toppling processes 
which occur in the model. When A = — 1 no transition 
is found between an active and an absorbing state, and 
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a diffusive dynamics of active sites is found. In the other 
extreme A = 1 the model corresponds to a conservative 
BTW sandpile, and our results indicate a discontinuous 
transition. We believe that the observed variations in the 
critical exponents may be due to crossover effects in the 
two limiting cases of the model. 

It is necessary to extend our simulations to values of 
A which are closer to both limits, in order to find out if 
the exponents approach limiting values. In particular, it 
would be interesting to find out if mean-field exponents 
are found as A — > — 1. The exponent ration v\\/v±_ should 
also be estimated through additional simulations, as well 
as the ratio m = (Pa)/Pa> whose value at £ = £ c is also 
universal. We believe the phase transition in the model 
to be discontinuous only at A = — 1, but the data we have 



collected so far does not allow us to discard the possibility 
that a tricritical point exists for some value of A between 
and 1. 

We are presently trying to answer these questions, so 
that the critical behavior of this model, which does not 
belong to the DP universality class, may be better known. 
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